#seed0=sample(1:2^15, 1)
seed0=4571
set.seed(seed0)



T=24
r=2
T0=240

M=200
N=1000

# tildeFM0[m,]=t(matrix(tildeF0,(T+T0)*r,1)) #t(tildeF[,1])
# rot1FM0[m,]=t(matrix(F0%*%H0,(T+T0)*r,1)) #t(tildeF[,1])
# rot2FM[m,]=t(matrix(F0%*%tildeH0,(T+T0)*r,1)) #t(tildeF[,1])


tildebarFM=matrix(NA,M,r)
rot1barFM=matrix(NA,M,r)
rot2barFM=matrix(NA,M,r)


tildeFM=matrix(NA,M,T*r)#=t(matrix(tildeF,T*r,1)) #t(tildeF[,1])
rot1FM=matrix(NA,M,T*r)#t(matrix(F%*%H,T*r,1)) #t(tildeF[,1])
rot2FM=matrix(NA,M,T*r)#t(matrix(F%*%tildeH,T*r,1)) #t(tildeF[,1])

tildeHM=matrix(NA,M,r*r)#=t(matrix(tildeF,T*r,1)) #t(tildeF[,1])
HM=matrix(NA,M,r*r)#t(matrix(F%*%H,T*r,1)) #t(tildeF[,1])

tildeLam=matrix(NA,N,r)
SigmaLam=matrix(NA,r,r)
SigmaLamM=matrix(NA,M,r)
SigmaLamtildeM=matrix(NA,M,r)


commtildeM=matrix(NA,M,(N*T))
commM=matrix(NA,M,(N*T))


#  when using T0

tildeLam0=matrix(NA,N,r)
SigmaLam0=matrix(NA,r,r)
SigmaLam0M=matrix(NA,M,r)
SigmaLam0tildeM=matrix(NA,M,r)





tildeFM0=matrix(NA,M,(T+T0)*r)#=t(matrix(tildeF,T*r,1)) #t(tildeF[,1])
rot1FM0=matrix(NA,M,(T+T0)*r)#t(matrix(F%*%H,T*r,1)) #t(tildeF[,1])
rot2FM0=matrix(NA,M,(T+T0)*r)#t(matrix(F%*%tildeH,T*r,1)) #t(tildeF[,1])

tildeHM0=matrix(NA,M,r*r)#=t(matrix(tildeF,T*r,1)) #t(tildeF[,1])
HM0=matrix(NA,M,r*r)#t(matrix(F%*%H,T*r,1)) #t(tildeF[,1])

tildebarFM0=matrix(NA,M,r)
rot1barFM0=matrix(NA,M,r)
rot2barFM0=matrix(NA,M,r)

tildebarFM0loc=matrix(NA,M,r)
rot1barFM0loc=matrix(NA,M,r)
rot2barFM0loc=matrix(NA,M,r)


tildeHM2=matrix(NA,M,r*r)#=t(matrix(tildeH2,r*r,1)) #t(tildeF[,1])
# HM2[m,]=t(matrix(H2,r*r,1)) #t(tildeF[,1])

tildeLam2M=matrix(NA,M,N*r)#   [m,]=t(matrix(tildeLam2,N*r,1)) #t(tildeF[,1])
#rot1FM[m,]=t(matrix(F%*%H,N*r,1)) #t(tildeF[,1])
rot2Lam2M=matrix(NA,M,N*r)#[m,]=t(matrix(Lambda%*%tildeH2,N*r,1)) #t(tildeF[,1])
rot2F2M=matrix(NA,M,T*r)#[m,]=t(matrix(F%*%t(solve(tildeH2)),T*r,1)) #t(tildeF[,1])
tildeF2M=matrix(NA,M,T*r)#[m,]=t(matrix(tildeF2,T*r,1))


tildeHM02=matrix(NA,M,r*r)#[m,]=t(matrix(tildeH02,r*r,1)) #t(tildeF[,1])
# HM2[m,]=t(matrix(H2,r*r,1)) #t(tildeF[,1])

tildeLam2M0=matrix(NA,M,N*r)#[m,]=t(matrix(tildeLam02,N*r,1)) #t(tildeF[,1])
#rot1FM[m,]=t(matrix(F%*%H,N*r,1)) #t(tildeF[,1])
rot2Lam2M0=matrix(NA,M,N*r)#[m,]=t(matrix(Lambda%*%tildeH02,N*r,1)) #t(tildeF[,1])
rot2F2M0=matrix(NA,M,T*r)#[m,]=t(matrix(F%*%t(solve(tildeH02)),T*r,1)) #t(tildeF[,1])
tildeF2M0=matrix(NA,M,T*r)#=t(matrix(tildeF02,T*r,1))




library(Matrix)
library(expm)

F=rnorm(T*r,0,1)
F=matrix(F,T,r)
F=F%*%solve(sqrtm(t(F)%*%F/T))

F00=rnorm(T0*r,0,1)
F00=matrix(F00,T0,r)
F00=F00%*%solve(sqrtm(t(F00)%*%F00/T0))


#SF=diag(r)
#SF[r,r]=r
SF=diag(abs(rnorm(r,1,1)))
F=F%*%SF



#V=(t(F)%*%F/T)
#SF%*%V%*%SF


#V[1:3,1:3]



m=1
for (m in 1:M)
{
  print(m)
  # estimation in small window T
  
  e=matrix(rnorm(T*N,0,1),T,N)
  Lambda=matrix(rnorm(r*N,0,1),N,r)
  X=F%*%t(Lambda)+e
  
  A=X%*%t(X)/(N*T)
  rr <- eigen(A)
  tildeF <- rr$vectors
  tildeF=sqrt(T)*tildeF[,1:r]
  mu_tildeF=apply(tildeF,2,mean)
  #tildeF=tildeF%*%diag(sign(mu_tildeF))
  lam <- rr$values
  lam=diag(lam) 
  R <- eigen(X%*%t(X)/(N*T)-diag(T)/T)
  Rlam <- R$values
  #R$vectors
  Uc=diag(sort(Rlam,decreasing=TRUE))
  Uc=Uc[1:r,1:r]
  tildeLam=t(X)%*%tildeF%*%solve(t(tildeF)%*%tildeF)
  
  
  
  
  
  
  tildeH=((t(Lambda)%*%Lambda)/N)%*%((t(F)%*%tildeF)/T)%*%solve(Uc)
  
  #tildeH
  #F%*%tildeH
  #tildeF
  
  
  
  SigmaF=(t(F)%*%F)/T
  rrr <- eigen(SigmaF)
  Ups <- rrr$vectors
  U <- rrr$values
  U=diag(U) 
  H=Ups%*%solve(sqrt(U))
  Q=sqrt(U)%*%t(Ups)
  
  #H=solve(Q)# same as Ups%*%solve(sqrt(U))
  #tildeH
  #H
  
  #F%*%H
  #F%*%tildeH
  #tildeF
  
  SigmaLam=solve(H)%*%(t(tildeLam)%*%tildeLam/N)%*%solve(t(H))
  SigmaLamM[m,]=diag(SigmaLam)
  
  
  SigmaLamtilde=solve(tildeH)%*%(t(tildeLam)%*%tildeLam/N)%*%solve(t(tildeH))
  SigmaLamtildeM[m,]=diag(SigmaLamtilde)
  
  
  tildeHM[m,]=t(matrix(tildeH,r*r,1)) #t(tildeF[,1])
  HM[m,]=t(matrix(H,r*r,1)) #t(tildeF[,1])
  
  tildeFM[m,]=t(matrix(tildeF,T*r,1)) #t(tildeF[,1])
  rot1FM[m,]=t(matrix(F%*%H,T*r,1)) #t(tildeF[,1])
  rot2FM[m,]=t(matrix(F%*%tildeH,T*r,1)) #t(tildeF[,1])
  
  tildebarFM[m,]=apply(tildeF,2,mean)
  rot1barFM[m,]=apply(F%*%H,2,mean)
  rot2barFM[m,]=apply(F%*%tildeH,2,mean)
  
  
  commtildeM[m,]=t(matrix(tildeF%*%t(tildeLam),T*N,1))
  commM[m,]=t(matrix(F%*%t(Lambda),T*N,1))
  
  #  estimation T but from loadings
  
  A2=t(X)%*%(X)/(N*T)
  rr2 <- eigen(A2)
  tildeLam2 <- rr2$vectors
  tildeLam2=sqrt(N)*tildeLam2[,1:r]
  mu_tildeLam=apply(tildeLam2,2,mean)
  tildeLam2=tildeLam2%*%diag(sign(mu_tildeLam))
  lam2 <- rr2$values
  lam2=diag(lam2) 
  R2 <- eigen(t(X)%*%(X)/(N*T)-diag(N)/N)
  Rlam2 <- R2$values
  #R$vectors
  Uc2=diag(sort(Rlam2,decreasing=TRUE))
  Uc2=Uc2[1:r,1:r]
  tildeF2=(X)%*%tildeLam2%*%solve(t(tildeLam2)%*%tildeLam2)
  
  
  
  
  
  
  tildeH2=((t(F)%*%F)/T)%*%((t(Lambda)%*%tildeLam2)/N)%*%solve(Uc2)
  
  # F times solve(t(tildeH2))
  # AF=F%*%t(solve(tildeH2))
  
  #tildeH
  #F%*%tildeH
  #tildeF
  
  
  
  # SigmaF=(t(F)%*%F)/T
  # rrr <- eigen(SigmaF)
  # Ups <- rrr$vectors
  # U <- rrr$values
  # U=diag(U) 
  # H=Ups%*%solve(sqrt(U))
  # Q=sqrt(U)%*%t(Ups)
  
  #H=solve(Q)# same as Ups%*%solve(sqrt(U))
  #tildeH
  #H
  
  #F%*%H
  #F%*%tildeH
  #tildeF
  
  #SigmaLam=solve(H)%*%(t(tildeLam)%*%tildeLam/N)%*%solve(t(H))
  #SigmaLamM[m,]=diag(SigmaLam)
  
  tildeHM2[m,]=t(matrix(tildeH2,r*r,1)) #t(tildeF[,1])
  # HM2[m,]=t(matrix(H2,r*r,1)) #t(tildeF[,1])
  
  tildeLam2M[m,]=t(matrix(tildeLam2,N*r,1)) #t(tildeF[,1])
  #rot1FM[m,]=t(matrix(F%*%H,N*r,1)) #t(tildeF[,1])
  rot2Lam2M[m,]=t(matrix(Lambda%*%tildeH2,N*r,1)) #t(tildeF[,1])
  rot2F2M[m,]=t(matrix(F%*%t(solve(tildeH2)),T*r,1)) #t(tildeF[,1])
  tildeF2M[m,]=t(matrix(tildeF2,T*r,1))
  # tildebarFM[m,]=apply(tildeF,2,mean)
  # rot1barFM[m,]=apply(F%*%H,2,mean)
  # rot2barFM[m,]=apply(F%*%tildeH,2,mean)
  
  
  
  # comment for now QQQ
  # estimation in large T0


  # DGP1 with lam0=0
  #  F0=rbind(F,matrix(0,T0,r))
  #  e0=matrix(rnorm(T0*N,0,1),T0,N)
  # # # #Lambda0=matrix(rnorm(r*N,0,1),N,r)
  #  X0=e0 #+F%*%t(Lambda)
  #  X0=rbind(X,X0)

  #
  # # DGP2 with lam0=0
  # F0=rbind(F,F00)
  # e0=matrix(rnorm(T0*N,0,1),T0,N)
  # Lambda0=matrix(rnorm(r*N,0,2),N,r)
  # X0=e0 +F00%*%t(Lambda0)
  # X0=rbind(X,X0)

  # # DGP3 with lam0=0
  F0=rbind(F,F00)
  e0=matrix(rnorm(T0*N,0,1),T0,N)
  Lambda0=matrix(rnorm(r*N,0,2),N,r)
  Lambda0[,r]=Lambda0[,r]/sqrt(N)
  X0=e0 +F00%*%t(Lambda0)
  X0=rbind(X,X0)



  A0=X0%*%t(X0)/(N*(T+T0))
  rr0 <- eigen(A0)
  tildeF0 <- rr0$vectors
  tildeF0=sqrt(T+T0)*tildeF0[,1:r]
  mu_tildeF0=apply(tildeF0,2,mean)
  tildeF0=tildeF0%*%diag(sign(mu_tildeF0))
  lam0 <- rr0$values
  lam0=diag(lam0)
  R0 <- eigen(X0%*%t(X0)/(N*(T+T0))-diag(T+T0)/(T+T0))
  Rlam0 <- R0$values
  #R0$vectors
  Uc0=diag(sort(Rlam0,decreasing=TRUE))
  Uc0=Uc0[1:r,1:r]

  SigmaF0=(t(F0)%*%F0)/(T+T0)
  rrr0 <- eigen(SigmaF0)
  Ups0 <- rrr0$vectors
  U0 <- rrr0$values
  U0=diag(U0)
  H0=Ups0%*%solve(sqrt(U0))
  Q0=sqrt(U0)%*%t(Ups0)



  tildeLam0=t(X0)%*%tildeF0%*%solve(t(tildeF0)%*%tildeF0)
  SigmaLam0=solve(H0)%*%(t(tildeLam0)%*%tildeLam0/N)%*%solve(t(H0))
  SigmaLam0M[m,]=diag(SigmaLam0)



  tildeH0=((t(Lambda)%*%Lambda)/N)%*%((t(F0)%*%tildeF0)/(T+T0))%*%solve(Uc0)

  SigmaLam0tilde=solve(tildeH0)%*%(t(tildeLam0)%*%tildeLam0/N)%*%solve(t(tildeH0))
  SigmaLam0tildeM[m,]=diag(SigmaLam0tilde)


  #tildeH
  #F%*%tildeH
  #tildeF




  #H=solve(Q)# same as Ups%*%solve(sqrt(U))
  #tildeH
  #H

  #F%*%H
  #F%*%tildeH
  #tildeF

  tildeHM0[m,]=t(matrix(tildeH0,r*r,1)) #t(tildeF[,1])
  HM0[m,]=t(matrix(H0,r*r,1)) #t(tildeF[,1])

  tildeFM0[m,]=t(matrix(tildeF0,(T+T0)*r,1)) #t(tildeF[,1])
  rot1FM0[m,]=t(matrix(F0%*%H0,(T+T0)*r,1)) #t(tildeF[,1])
  rot2FM0[m,]=t(matrix(F0%*%tildeH0,(T+T0)*r,1)) #t(tildeF[,1])

  tildebarFM0[m,]=apply(tildeF0,2,mean)
  rot1barFM0[m,]=apply(F0%*%H0,2,mean)
  rot2barFM0[m,]=apply(F0%*%tildeH0,2,mean)

  tildebarFM0loc[m,]=apply(tildeF0[1:T,],2,mean)
  rot1barFM0loc[m,]=apply(F0[1:T,]%*%H0,2,mean)
  rot2barFM0loc[m,]=apply(F0[1:T,]%*%tildeH0,2,mean)

  # > cor(tildeFM0[,1],rot2FM0[,1])
  # [1] 0.9687704
  # > cor(tildeFM[,1],rot2FM[,1])
  # [1] 0.9686775
  # > cor(tildeFM0[,24],rot2FM0[,24])
  # [1] 0.9952786
  # > cor(tildeFM[,24],rot2FM[,24])
  # [1] 0.9952997



  #    cor(tildeFM0[,24],rot2FM0[,24])
  #    cor(tildeFM[,25],rot2FM[,25])
  #
  #  from loadings but T0




  A02=t(X0)%*%(X0)/(N*(T+T0))
  rr02 <- eigen(A02)
  tildeLam02 <- rr02$vectors
  tildeLam02=sqrt(N)*tildeLam02[,1:r]
  mu_tildeLam02=apply(tildeLam02,2,mean)
  tildeLam02=tildeLam02%*%diag(sign(mu_tildeLam02))
  lam02 <- rr02$values
  lam02=diag(lam02)
  R02 <- eigen(t(X0)%*%(X0)/(N*(T+T0))-diag(N)/(N))
  Rlam02 <- R02$values
  #R0$vectors
  Uc02=diag(sort(Rlam02,decreasing=TRUE))
  Uc02=Uc02[1:r,1:r]

  # SigmaF0=(t(F0)%*%F0)/(T+T0)
  # rrr0 <- eigen(SigmaF0)
  # Ups0 <- rrr0$vectors
  # U0 <- rrr0$values
  # U0=diag(U0)
  # H0=Ups0%*%solve(sqrt(U0))
  # Q0=sqrt(U0)%*%t(Ups0)
  #


  #   tildeLam0=t(X0)%*%tildeF0%*%solve(t(tildeF0)%*%tildeF0)
  #    SigmaLam0=solve(H0)%*%(t(tildeLam0)%*%tildeLam0/N)%*%solve(t(H0))
  #    SigmaLam0M[m,]=diag(SigmaLam0)



  tildeH02=((t(F0)%*%F0)/(T+T0))%*%((t(Lambda)%*%tildeLam02)/(N))%*%solve(Uc02)
  #tildeH
  #F%*%tildeH
  #tildeF
  tildeF02=(X0)%*%tildeLam02%*%solve(t(tildeLam02)%*%tildeLam02)






  # tildeH2=((t(F)%*%F)/T)%*%((t(Lambda)%*%tildeLam2)/N)%*%solve(Uc2)

  # F times solve(t(tildeH2))
  # AF=F%*%t(solve(tildeH2))

  #tildeH
  #F%*%tildeH
  #tildeF



  # SigmaF=(t(F)%*%F)/T
  # rrr <- eigen(SigmaF)
  # Ups <- rrr$vectors
  # U <- rrr$values
  # U=diag(U)
  # H=Ups%*%solve(sqrt(U))
  # Q=sqrt(U)%*%t(Ups)

  #H=solve(Q)# same as Ups%*%solve(sqrt(U))
  #tildeH
  #H

  #F%*%H
  #F%*%tildeH
  #tildeF

  #SigmaLam=solve(H)%*%(t(tildeLam)%*%tildeLam/N)%*%solve(t(H))
  #SigmaLamM[m,]=diag(SigmaLam)

  tildeHM02[m,]=t(matrix(tildeH02,r*r,1)) #t(tildeF[,1])
  # HM2[m,]=t(matrix(H2,r*r,1)) #t(tildeF[,1])

  tildeLam2M0[m,]=t(matrix(tildeLam02,N*r,1)) #t(tildeF[,1])
  #rot1FM[m,]=t(matrix(F%*%H,N*r,1)) #t(tildeF[,1])
  rot2Lam2M0[m,]=t(matrix(Lambda%*%tildeH02,N*r,1)) #t(tildeF[,1])
  rot2F2M0[m,]=t(matrix(F%*%t(solve(tildeH02)),T*r,1)) #t(tildeF[,1])
  tildeF2M0[m,]=t(matrix(tildeF02[1:T,],T*r,1))




  #H=solve(Q)# same as Ups%*%solve(sqrt(U))
  #tildeH
  #H

  #F%*%H
  #F%*%tildeH
  #tildeF

  # tildeHM0[m,]=t(matrix(tildeH0,r*r,1)) #t(tildeF[,1])
  # HM0[m,]=t(matrix(H0,r*r,1)) #t(tildeF[,1])
  #
  # tildeFM0[m,]=t(matrix(tildeF0,(T+T0)*r,1)) #t(tildeF[,1])
  # rot1FM0[m,]=t(matrix(F0%*%H0,(T+T0)*r,1)) #t(tildeF[,1])
  # rot2FM0[m,]=t(matrix(F0%*%tildeH0,(T+T0)*r,1)) #t(tildeF[,1])
  #
  # tildebarFM0[m,]=apply(tildeF0,2,mean)
  # rot1barFM0[m,]=apply(F0%*%H0,2,mean)
  # rot2barFM0[m,]=apply(F0%*%tildeH0,2,mean)
  #
  # tildebarFM0loc[m,]=apply(tildeF0[1:T,],2,mean)
  # rot1barFM0loc[m,]=apply(F0[1:T,]%*%H0,2,mean)
  # rot2barFM0loc[m,]=apply(F0[1:T,]%*%tildeH0,2,mean)
  
  
}


mm=1
rot1FMM=matrix(NA,M,T*r)
rot1FMM0=matrix(NA,M,(T+T0)*r)

for (mm in  (1:M))
{
  
  rot1FMM[mm,]=abs(rot1FM[mm,])*sign(apply(tildeFM,2,mean))
  rot1FMM0[mm,]=abs(rot1FM0[mm,])*sign(apply(tildeFM0,2,mean))
  
  
}

apply(SigmaLamM,2,mean)/apply(SigmaLam0M,2,mean)

apply(SigmaLamtildeM,2,mean)/apply(SigmaLam0tildeM,2,mean)




apply(tildebarFM,2,mean)
rot1barFM[1,]
100*(apply(rot2barFM,2,mean)/apply(tildebarFM,2,mean) -1 )



apply(tildebarFM0,2,mean)
rot1barFM0[1,]
apply(rot2barFM0,2,mean)
100*(apply(rot2barFM0,2,mean)/apply(tildebarFM0,2,mean) -1 )


c1=0
c10=0
c2=0
c20=0

m=1
for (m in 1:T)
{
  c10=c10+cor(tildeFM0[,m],rot2FM0[,m])
  # [1] 0.9687704
  c1=c1+cor(tildeFM[,m],rot2FM[,m])
  c20=c20+cor(tildeFM0[,T+T0+m],rot2FM0[,T+T0+m])
  # [1] 0.9687704
  c2=c2+cor(tildeFM[,T+m],rot2FM[,T+m])
  
  print(cbind(cor(tildeFM0[,m],rot2FM0[,m]),cor(tildeFM[,m],rot2FM[,m]),cor(tildeFM0[,T+T0+m],rot2FM0[,T+T0+m]),cor(tildeFM[,T+m],rot2FM[,T+m]),cor(commM[,m],commtildeM[,T+m]))
        
        
        )
  
  # [1] 0.9686775
  #cor(tildeFM0[,T],rot2FM0[,T])
  # [1] 0.9952786
  #cor(tildeFM[,T],rot2FM[,T])
  # [1] 0.9952997
  #cor(tildeFM0[,T+1],rot2FM0[,T+1])
  # [1] 0.9952786
  #cor(tildeFM[,T+1],rot2FM[,T+1])
  
  #cor(tildeFM0[,T+T0],rot2FM0[,T+T0])
  
}
c1/T
c10/T
c2/T
c20/T


c10full=0
c20full=0
m=1
for (m in 1:(T+T0))
{
  c10full=c10full+cor(tildeFM0[,m],rot2FM0[,m])
  c20full=c20full+cor(tildeFM0[,T+T0+m],rot2FM0[,T+T0+m])
}
c10full/(T+T0)
c20full/(T+T0)



hist(tildeFM[,1],breaks=100)
abline(v=rot1FMM[1,1],col=2)
abline(v=mean(tildeFM[,1]),col=3)
hist((tildeFM[,1]- mean(tildeFM[,1]))/sd(tildeFM[,1]),breaks=100,freq=FALSE,prob=TRUE)
curve(dnorm(x, mean=0, sd=1),add=TRUE)

pdf('tildeFM11.pdf')
tildeFM11 <- hist((tildeFM[,1]- mean(tildeFM[,1]))/sd(tildeFM[,1]),breaks=100, plot=FALSE) #hist(mfunds$tbill , breaks = 30 , plot = FALSE)
plot(tildeFM11,  main="HISTOGRAM: LOCAL PCA (T=2, r=2)", xlab="F(1,1)",freq=FALSE) # Plot 1st histogram using a transparent color
curve(dnorm(x, mean=0, sd=1), lwd=2, add=TRUE)
dev.off()

pdf('tildeFM12.pdf')
tildeFM12 <- hist((tildeFM[,2]- mean(tildeFM[,2]))/sd(tildeFM[,2]),breaks=100, plot=FALSE) #hist(mfunds$tbill , breaks = 30 , plot = FALSE)
plot(tildeFM12,  main="HISTOGRAM: LOCAL PCA (T=2, r=2)", xlab="F(1,2)",freq=FALSE) # Plot 1st histogram using a transparent color
curve(dnorm(x, mean=0, sd=1), lwd=2, add=TRUE)
dev.off()

pdf('tildeFM21.pdf')
tildeFM21 <- hist((tildeFM[,3]- mean(tildeFM[,3]))/sd(tildeFM[,3]),breaks=100, plot=FALSE) #hist(mfunds$tbill , breaks = 30 , plot = FALSE)
plot(tildeFM21,  main="HISTOGRAM: LOCAL PCA (T=2, r=2)", xlab="F(2,1)",freq=FALSE) # Plot 1st histogram using a transparent color
curve(dnorm(x, mean=0, sd=1), lwd=2, add=TRUE)
dev.off()

pdf('tildeFM22.pdf')
tildeFM22 <- hist((tildeFM[,4]- mean(tildeFM[,4]))/sd(tildeFM[,4]),breaks=100, plot=FALSE) #hist(mfunds$tbill , breaks = 30 , plot = FALSE)
plot(tildeFM22,  main="HISTOGRAM: LOCAL PCA (T=2, r=2)", xlab="F(2,2)",freq=FALSE) # Plot 1st histogram using a transparent color
curve(dnorm(x, mean=0, sd=1), lwd=2, add=TRUE)
dev.off()

hist(tildeFM[,2],breaks=100)
abline(v=rot1FMM[1,2],col=2)
abline(v=mean(tildeFM[,2]),col=3)


hist((tildeFM[,2]- mean(tildeFM[,2]))/sd(tildeFM[,2]),breaks=100,prob=TRUE)
curve(dnorm(x, mean=0, sd=1),add=TRUE)



hist(tildeFM[,T*r],breaks=100)
abline(v=rot1FMM[1,T*r],col=2)
abline(v=mean(tildeFM[,T*r]),col=3)
hist((tildeFM[,T*r]- mean(tildeFM[,T*r]))/sd(tildeFM[,T*r]),breaks=25,prob=T)
curve(dnorm(x, mean=0, sd=1),add=TRUE)



hist(tildeFM[,3],breaks=100)
abline(v=rot1FMM[1,3],col=2)
abline(v=mean(tildeFM[,3]),col=3)

hist((tildeFM[,3]- mean(tildeFM[,3]))/sd(tildeFM[,3]),breaks=25,prob=T)
curve(dnorm(x, mean=0, sd=1),add=TRUE)

